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Abstract 

The purpose of this paper is to provide analytical and numerical solutions of the forma¬ 
tion and evolution of the localized plastic zone in a uniaxially loaded bar with variable 
cross-sectional area. An energy-based variational approach is employed and the governing 
equations with appropriate physical boundary conditions, jump conditions, and regular¬ 
ity conditions at evolving elasto-plastic interface are derived for a fourth-order explicit 
gradient plasticity model with linear isotropic softening. Four examples that differ by 
regularity of the yield stress and stress distributions are presented. Results for the load 
level, size of the plastic zone, distribution of plastic strain and its spatial derivatives, plas¬ 
tic elongation, and energy balance are constructed and compared to another, previously 
discussed non-variational gradient formulation. 
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1. Introduction 

The presence of a softening branch of the stress-strain curve, usually caused by ini¬ 
tiation, propagation and coalescence of defects such as micro-cracks or micro-voids, is 
a phenomenon typical of quasi-brittle materials. Softening often leads to localization of 
strain into narrow bands whose width is related to an intrinsic length dictated by the het¬ 
erogeneities of the material microstructure. Softening can be conveniently incorporated 
into damage models, but can also be described by plasticity with a negative hardening 
modulus. However, constitutive models within the classical continuum framework of sim¬ 
ple materials do not contain any length scale reflecting the typical size of microstructural 
features. Therefore, the localization processes due to softening are not described prop¬ 
erly and mathematical models lead to ill-posed problems accompanied by localization of 
strain into subdomains of zero volume and consequently to vanishing dissipation. Various 
enrichments incorporating some information about the material heterogeneity have been 
developed. They utilize, for example, additional kinematic variables, weighted spatial av¬ 
erages, higher-order gradients, or rate-dependent terms; see e.g. the comparative studies 
and review papers by de Vree et al. (1995), Jirasek (1998), Peerlings et al. (2001), Jirasek 
and Rolshoven (2003), Bazant and Jirasek (2002), Jirasek and Rolshoven (2009a) and 
Jirasek and Rolshoven (2009b). These techniques on the one hand preclude localization 
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and loss of ellipticity of the governing equations, but on the other hand significantly com¬ 
plicate the overall analysis. For instance, in the case of higher-order gradient models, the 
regularity conditions of internal variables at the evolving elasto-plastic interface are not 
easy to characterize. 

The present paper is devoted to the investigation of the Aifantis explicit fourth-order 
gradient plasticity model, cf. e.g. Zbib and Aifantis (1988) or Miihlhaus and Aifantis 
(1991), under conditions leading to non-uniform stress fields. To keep the analysis trans¬ 
parent, we confine ourselves to one-dimensional tension tests and perform our analysis 
in the framework of the so-called energetic solutions introduced in an abstract setting 
in Mielke and Theil (2004) and in the context of finite-strain plasticity in Mielke (2003), 
generalizing earlier variational formulations of damage (Francfort and Marigo, 1993) and 
fracture (Francfort and Marigo, 1998). Building on this basis, we will derive the govern¬ 
ing equations and appropriate boundary and jump conditions. In particular, the often 
questioned regularity conditions for internal variables at the elasto-plastic interface will 
emerge naturally in a consistent and unified way. Detailed solutions for four different 
test problems with various regularity of the yield stress and stress distributions will be 
presented and compared to results obtained from the standard non-variational gradient 
formulation available in Jirasek et al. (2010). The present work can also be viewed as a 
continuation and extension of results provided in our previous paper (Jirasek et ah, 2013), 
where the second-order gradient plasticity model was investigated. The main difference 
from our previous work is that now we treat a more complicated model with higher-order 
regularity requirements on internal variables and more intricate conditions at the elasto- 
plastic interface. In addition, we approach the problem utilizing the nowadays standard 
variational framework for rate-independent evolution. 

This study is also closely related to several one-dimensional studies into energy-based 
second-order gradient models of strain-softening damage and plasticity. In particular, 
Pham et ah (2011) performed a detailed analysis of stability and bifurcation of localized 
and homogenous states, obtained in the closed form, for a parameterized family of gradient 
damage models. These results were later refined by Pham and Marigo (2013), who studied 
size effects and snap-back behaviour at the structural scale predicted by the same group of 
damage models. The combination of damage and plasticity has been the subject of recent 
contributions by Del Piero et ah (2013) and Alessi et ah (2014), with the emphasis on the 
competition between brittle and ductile failure; the former work, as well as Milasinovic 
(2004), also contain a validation against experimental data. Theoretical results of fracture 
and plasticity as T-limits of damage models within one-dimensional setting are described 
in Iurlano (2013), and non-local damage or fracture analyses in bars under tension are 
discussed in Jirasek and Zeman (2015) and Lellis and Royer-Carfagni (2001); all these 
works also employ a variational formulation. Our analysis extends these contributions 
by treating a model regularized by the fourth derivative of internal variables and by 
deriving the regularity of internal variables directly from energy-minimization arguments, 
rather than enforcing them through additional boundary conditions at the moving elastic- 
inelastic interfaces. 

The energetic formulation for rate-independent processes comprises several steps and 
relies on two principles. In the abstract setting (Mielke, 2006), the state of the system 
within a fixed time horizon T is described in terms of a ”non-dissipative” field u(t,x), 
u(t ) G U, x G D, where D denotes the spatial domain, and a ’’dissipative” field z(t,x), 
z(t ) G Z, which specifies the irreversible processes at time t G [0, T\. The state of 
the system is fully characterized by the state variables q(t,x), q(t ) = ( u(t),z(t )) G Q = 


2 


u x Z. Typically, u is the displacement field and 2 is the field of internal variables related 
to the inelastic phenomena, such as plastic strain or damage. Further, we consider the 
total free (Helmholtz) energy of the body £ : [ 0 , T] x Q —> M U {+00} together with the 
dissipation distance T> : Z x Z —>■ M U {+00} which specifies the minimum amount of 
energy spent by the continuous transition from state 2 A to state 2 ( - 2 \ For notational 
convenience, we will sometimes refer to the dissipation distance by T>{qd\ q^ 2 ’) instead of 
'D{z <yl \ 2 ^). Then, the process q : [0, T] —> Q is an energetic solution to the initial-value 
problem described by (£,2?,q 0 ) if it satisfies 

(i) Global stability: for all t E [ 0 , T] and for all q G Q 

£(t,q(t)) < £(t,q)+ V(q(t),q) (S) 


which ensures that the solution minimizes the sum £ + V, 

(ii) Energy equality: for all t G [0, T] 

£(t,q(t)) + VarD(g;0,t) = £(0,g(0)) + f V(s)ds (E) 

Jo 

which expresses energy balance in terms of the internal energy, dissipated energy Varx>, 
and time-integrated power of external forces V, 

(iii) Initial condition: 

9 ( 0 ) = <lo (!) 

The dissipation along a process q is expressed as 


Var£>(qr; 0, t ) = sup 




2—1 


( 1 ) 


where the supremum is taken over all neN and all partitions of the time interval [0, t], 
0 — to < ti <■■■< t n — t. Together, the two principles (S) and (E) along with initial 
condition (I) naturally give rise to an 

Incremental problem: for k = 1,..., N 


q(tk ) e Arg min [£(t k ,q) + V(q(t k -i),q)] (IP) 

5eQ 

amenable to a numerical solution, in which each step is realized as a minimization prob¬ 
lem, e.g. Ortiz and Stainier (1999); Carstensen et al. (2002); Petryk (2003). The main 
conceptual difficulty with this incremental problem is that it represents a global min¬ 
imization, which is computationally cumbersome and physically difficult to justify for 
non-convex energies. It is reasonable, however, to assume that stable solutions to (IP) 
are associated with local minima; for comparative studies into evolution driven by local 
and global energy minimization see e.g. Mielke (2011); Braides (2014); Roubfcek (2015). 
On the other hand, the variational approach offers many advantages, among which we 
highlight that it provides a unified setting for the analysis, allows for discontinuities in 
space, incorporates the governing laws with boundary conditions and provides regularity 
conditions at the elasto-plastic interface. 

For the uniaxial displacement-controlled tension test, we can further specify all the 
quantities introduced above in more detail. Displacement u(t,x), where we have used 
the light face letter since u is now a scalar field, as a function of the spatial coordi¬ 
nate 1 6 O C K, represents the ”non-dissipative” component; the total linearized strain is 
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simply the spatial derivative u'(t, x ) = du(t, x)/dx. The ’’dissipative” variables, describing 
the irreversible processes, are plastic strain e p (t,x) and cumulative plastic strain n{t,x), 
i.e. z(t,x) = (e p (t, x), n(t, a;)). 1 The corresponding function spaces are as follows: 

u E 14(f) = {«6 W 1,2 (Q) | u = unit) on dVL in the sense of traces} (2a) 

% g v ep = w 2 ’ 2 (n ) ( 2 b) 

KEV K = {keW 2 ' 2 (Q)\k>0} (2c) 

where Uu(t) denotes prescribed displacements on the boundary (specifying the Dirichlet 
boundary condition) and W k ’ 2 stands for the space of all Lebesgue square-integrable 
functions with square-integrable generalized derivatives up to order k\ later on, we will 
employ a subset Wq ’ 2 consisting of functions vanishing at the boundary, for details we 
refer to e.g. Evans (2010). Consequently, we identify U = V u (f), which now depends 
on time (due to the time-dependent values presented on the boundary), Z = V £ x V K , 
Q = U x Z = 14(f) x V e x V K , and define the total free energy of the body 

£(t,u,e p , K,) = I ~EA(u'— E p ) 2 dx + f -HA{k 2 — I a k" 2 ) dx — f Abudx (3) 
Jn 2 Jn 2 Jo 

and the dissipation distance 


V(z^\z^) 



+ oo otherwise 


Quantities appearing in the definitions of energies represent the Young modulus E [Pa], 
softening modulus H < 0 [Pa], characteristic length of the material l [m], initial yield 
stress <j 0 [Pa], a function describing the distribution of the cross-sectional area along the 
bar A(x) [m 2 ], and prescribed body force density b{x) [N/m 3 ]. For completeness, let us 
note that the power of external forces in (E) has the form V(s) = F(s)ud{s), where F(s) 
denotes the reaction force as a function of time and the dot stands for the time derivative. 

The paper is organized as follows. In Section 2, we will revisit the energetic formula¬ 
tion for the case of monotone loading, i.e. uo(t) A 0, which greatly simplifies the specific 
form of (S), (E), (I), and (IP) accompanied by (2)-(4). Further, the governing equations 
with boundary conditions, jump conditions, and regularity conditions at the elasto-plastic 
interface will be derived for the resulting one-dimensional fourth-order gradient-enriched 
plasticity model. Sections 3 and 4 are concerned with piecewise constant yield stress and 
stress distributions, which also represent the only cases amenable to analytical solutions. 
There, it will be shown that the structural response may, in a certain range, exhibit 
hardening due to the gradient enrichment despite the softening character of the material 
model. In Sections 5 and 6, the results for piecewise linear and quadratic stress held 
distributions will be compared to standard non-variational solutions available in Jirasek 
et al. (2010). In spite of all the simplifications we will be forced to use numerical so¬ 
lutions. The influence of data variation to the evolution of the plastic zone, its profile 
and load-displacement diagrams will also be investigated. Finally, in Appendix A, the 


1 Strictly speaking, the ”non-dissipative” component should have read u e i(x, t) = u(x, t) — f e p [x, t ) dx, 

where J • dx denotes a primitive integral of a function •. Since such an affine transformation does 
not affect the solution, we adopt u instead of u e i as our primal variable in further considerations for 

convenience. 
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stability conditions for the case of a uniform bar are discussed, and in Appendix B, we 
verify optimality of the obtained regularity conditions at the elasto-plastic interface by 
an independent argument. 


2. Energy-Based Formulation 

2.1. General Considerations 

For S and T> given by (3) and (4), minimization in (IP) with respect to k gives K(tk) = 
k(4_i) + | E p — £ p (t k _ i)|. Consequently, (IP) reduces to 


(u(t k ),£ p (t k )) £ Arg min / -EA{u-e p ) 2 Ax 


(u,£ p )eV u (tk)xV Sp Jn 


1 


+ / -HA[K(t k _!) + \e p - £p(4_i)|] 2 dx 
Jn * 

1 — ‘ / 4 r // 


( 5 ) 


in 


AI a {k"( tk-x) + I £ p - £ p (t k -i)\'} 2 Ax 


+ / A(Jq\e p — £ p (t k -i)\Ax — / AbuAx 
J £1 J 


cf. also Section 4.3 in Mielke (2003), where the local plasticity model with hardening is 
discussed. In Eq. (5), spatial derivatives are understood in the sense of distributions. For 
further considerations, we will restrict ourselves to tensile loading with possible clastic 
unloading, but never with a reversal of the plastic flow. Then, the plastic strain e p and 
the cumulative plastic strain k are equal, and we can use k as the only internal variable. 
As a result, instead of the incremental approach given in Eqs. (IP) and (5), it is fully 
sufficient to consider a total formulation providing a parameterized solution (■ u(t),K,(t )) 
which does not violate the irreversibility constraints. Note that such a parametrization 
is mathematically justified only when the elastic energy £ is strictly convex in q. which 
implies that the solution is time-continuous for sufficiently regular loading, cf. Mielke 
and Tlieil (2004). Later on, instead of imposing the Dirichlet boundary conditions, we 
prescribe directly the size of the plastic zone as a function of time in order to control 
the system evolution. This approach automatically entails time-continuity of q(t), which 
implies satisfaction of the energy balance (Pham et ah, 2011; Pham and Marigo, 2013), 
and also justifies the total formulation. 

Taking into account all the above simplifications, the minimization problem in Eq. (5) 
reduces to 

(u(t),n(t)) E Arg min II (u,k) (6) 

(u,K)&Vu(t)x \4 


where 


II(it, 7c) = / -EA{u — n) 2 Ax + 
.In ^ 


-HA(k 2 - l 4 ?i" 2 ) dx 


Ao-qk dx 


AbuAx 


( 7 ) 


Jn Jn 

which resembles a variational inequality of the first kind, due to the requirement k > 0 
in Eq. (2c). 
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The first variation (Gateaux derivative) of II furnishes us with the optimality condition 


<HI(u, k; 5u, 5 k) = / EA(u' — k)(5u'— 5 k) dx + / HA(k5k — in" 8 k”) dx 


'O 


>n 


( 8 ) 


+ / Axo<5k dx — Abdu dx > 0 

Jn Jn 


where 5u and 5k are admissible variations satisfying (u + 5u, k + 5k) € I4(t) x 14 (for 
simplicity, explicit dependence on time t has been dropped). Employing the integration 
by parts we arrive at 


5H(u, k]5u,5k) = — f \{EA{v! — k )) 1 + Ab]5u dx 

Jn 

+ [ [HAk - (j HAI a k .")" + Aa 0 - EA{u’ - «)]<5«da; 

Jn 

+ EA(y! — k)ti5u — \EA(u — K)5u\ Xi 
sn i 

- HAl A K"n5K' + YJHAI a k"5k%. 
sn i 

+ ^( HAl A K")'n5K - YJ{HAl i K")'5K\ Xi 
sn i 

where denotes the boundary integral, in our one-dimensional setting reduced to the 
sum over two end points of the interval fl, and n is the unit outer normal, which equals —1 
at the left and 1 at the right dQ r part of the boundary dQ = U dVt^. The sums 
YJi are taken over all points of possible discontinuity Xi and 

Iflxi = f(xt) ~ f(xj) = lim f( x ) - lim f(x) (10) 

X^Xi XJXi 

represents the jump of function f(x) at x t . The necessary condition for a local minimum 
is non-negativity of the first variation of the functional II for all admissible variations 5u 
and 5k, cf. Section 8.4.2 in Evans (2010), or Chapter 5 in Roubicek (2010). Below we show 
that such an approach leads to a consistent set of governing equations, namely the equi¬ 
librium equations, complementarity conditions of the plastic flow, boundary conditions, 
and regularity conditions at the elasto-plastic interfaces. Analysis of the second varia¬ 
tion (second-order Gateaux derivative), which is related to stability of the solution (S), is 
postponed to Appendix A. 

2.2. Governing Equations 

Since u + 5u G 14, we have 5u G W^' 2 {Gl), meaning that 5u is arbitrary inside fl with 
zero trace on the physical boundary dfl. Thus the expression multiplying 5u in the first 
line of (9) must vanish, providing us with the static equilibrium condition 

( EA(u' — k ))' + Ab = 0 in O (11) 

Here, v! corresponds to the total strain, u' — k is the elastic strain, and EA(u' — k) is 
the axial force, which is required to be continuous according to the third line of (9), 
because 5u is arbitrary inside fh Due to the zero trace of 5u, the sum over boundary 
points in the third line of (9) vanishes. 
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Since k + 8k G V K , variations 8k cannot be completely arbitrary. Let us define the 
plastic zone as the open set X p = {x E 1 k(x) > 0}, i.e., as the support of ft, and 
the elastic zone as the open set I e = {x G fi | k(x) = 0}. In X p , the variation 8k can 
have an arbitrary sign, and so the expression multiplying 8k in the second line of (9) 
must vanish. On the other hand, only non-negative variations 8k are admissible in X e , 
and so the expression multiplying 8k does not necessarily vanish but is constrained to be 


non-negative. The resulting conditions 

HAk - ( HAI a k ")" + A(j q = EA(u' - k) in X p ( 12 ) 

HAk - (. HAI a k ")" + Aa 0 > EA(u' - k) in X e (13) 

combined with the definitions of X p and X e can be presented in the complementarity format 

k > 0 (14) 

HAk - (. HA1 4 k ")" + Aa 0 - EA(v! - k) > 0 (15) 

[HAk - ( HA1 4 k ")" + Aa 0 - EA(v! - ft)] ■ ft = 0 (16) 

Note also that since k — 0 in X e , condition (13) could be simplified to 

Acr 0 > EAv! in X e (17) 


For AH = const., the second term on the left-hand side of Eq. (12) reduces to 
—LL4Z 4 ft IV , and the standard formulation of the fourth-order gradient plasticity model 
is recovered, cf. Jirasek et al. (2010) and Tab. 1. For variable sectional area and/or 
variable plastic modulus, expansion of the second term on the left-hand side of Eq. (12) 
gives — 1 4 [H(x)A(x)k w (x) + 2(H(x)A{x))' k"'(x) + (H (x)A(x))" k" (x)\. With increasing 
magnitude of the derivatives of H{x)A{x) we expect also increasing differences between 
the solutions corresponding to the classical and variational formulations. 

In addition to conditions (11)—(13), which have been deduced as optimality conditions 
following from the first two lines of (9), the last two lines of (9) provide us with boundary 
and regularity conditions for the plastic strain. 

Let us first discuss the boundary conditions. Again, we have to distinguish between 
the plastic part of the physical boundary, dQ P\X p , and the elastic part, dQ fl X e . 

1. Boundary point in a plastic state, characterized by k > 0: 

The variation 8k as well as its derivative 8k' at such a point can have an arbi¬ 
trary sign and the terms that multiply them must vanish. This leads to boundary 
conditions k" = 0 and (HAI a k")' = 0. 

2. Boundary point in an elastic state, characterized by k = 0: 

The variation 8k at such a point can only be zero or positive, and so the term that 
multiplies 8k in the first sum in the fifth line of (9) must not be negative but does 
not need to vanish. Therefore, for a boundary point in an elastic state we obtain 
the inequality condition (HAl A k")' n > 0. Recall that n = —1 at the left boundary 
and n — 1 at the right boundary. Regarding the second condition, tested by the 
derivative of the variation of plastic strain, we have to distinguish the following two 
subcases: 

(a) Nonzero derivative of plastic strain at the boundary: 

If ft = 0 and ft' 7 ^ 0 at a boundary point (which necessarily means ft' > 0 at the 
left boundary and ft' < 0 at the right boundary, by virtue of the universally 
valid admissibility condition k(x) > 0 ), then the variation of plastic strain, 
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5k', has an arbitrary sign, and the term HA1 4 k"u that multiplies 5 k' in the 
first sum in the fourth line (9) must vanish. Since HAl 4 ^ 0, this gives the 
boundary condition k" — 0 . 

(b) Zero derivative of plastic strain at the boundary: 

If k = 0 and k' = 0 at a boundary point, then the derivative of the variation 
of plastic strain, 5k', can be positive at the left boundary and negative at 
the right boundary, which means that Sk'ti can be negative. Consequently, 
the term HA1 4 k" that multiplies 5k' n in the first sum in the fourth line (9) 
must not be negative (note the minus sign before the sum). Since H < 0 and 
l 4 A > 0, the resulting condition reads k" < 0. But this actually cannot be 
satisfied as a strict inequality, because k" < 0 in combination with k = 0 and 
k' — 0 would lead to a violation of the admissibility condition k(x) > 0 in the 
near vicinity of the boundary point. So once again, we conclude that k" must 
vanish. 

We have found that the boundary condition k" = 0 applies independently of the state of 
the material at the boundary. On top of that, we have ( HA1 4 k ")' = 0 and k > 0 if the 
boundary point is in a plastic state, and ( HAl 4 K")'n > 0 and k = 0 if the boundary point 
is in an elastic state. All this can be summarized by the following boundary conditions: 

k" = 0 (18) 

(H Al 4 k")' n > 0 (19) 

k > 0 ( 20 ) 

(H Al 4 k")' nK = 0 (21) 

Let us now turn our attention to the continuity or regularity conditions that 
can be deduced from the jump terms in the last two lines of (9). Again, in the plastic 
domain the variation of plastic strain and its derivative have arbitrary signs and remain 
continuous, and so the jumps in HA1 4 k" and in ( HA1 4 k ")' must vanish. In other words, 
continuity of these terms must be preserved. Note that k and k' are continuous by 
assumption, but continuity of k" or k'" is not assumed apriori, and in fact is not maintained 
at points where for instance H or A have a jump. 

Inside the elastic domain, the plastic strain is identically zero and thus all its deriva¬ 
tives are zero, too, which means that the corresponding jump terms in the last two lines 
of (9) automatically vanish. However, special attention should be paid to those points 
of the boundary of the elastic domain which at the same time belong to the closure of 
the plastic domain (recall that the plastic domain is an open set), i.e., to the points of 
the elastoplastic interface, formally defined as dX ep = Z e D X p . Since these are not 
internal points of X e , we cannot directly infer that all derivatives of k vanish here. Con¬ 
tinuous differentiability of k implies that k = 0 and k! — 0 at dX ep , but higher derivatives 
could in principle exhibit a jump. So it is necessary to examine again the corresponding 
jump terms in (9). The variation 5k at dX ep can be zero or positive, but never negative. 
Therefore, the resulting optimality condition is the inequality \(HA1 4 k")'\ < 0. 

The last condition to be derived is the most delicate one. The derivative of the vari¬ 
ation of plastic strain, 5k', cannot be set to a nonzero value at a point of dX ep without 
simultaneously prescribing a positive value of 5k, otherwise the admissibility condition 
k(x) + 5k(x) > 0 would be violated in the vicinity of that point. Still, various combi¬ 
nations of 5k and 5k' can be selected such that the latter becomes increasingly “more 


Table 1: Comparison of standard and variational formulation: governing equations, 
boundary conditions, and regularity conditions for internal variable k. 


Standard formulation 

Variational formulation 

Note 

HAk - HA1 4 k w + 

+Aa 0 = EA(u' — k), k > 0 

HAk - ( HA1 4 k")"+ 

-\-Aoq = EA(u' — k), k > 0 

in X p 

Aa 0 > EAu', k = 0 

Aa 0 > EAu', k = 0 

in X e 

k > 0, k" = 0, k'" = 0 

o 

\ 

o' 

o' 

A 

£ 

at dfl D X p 

O 

o' 

£ 

k = 0 , k" = 0, ( HAl 4 K")'n > 0 

at on n X e 

continuous k, k!, k" 

continuous k, k!, HA1 4 k" 

in 

continuous k!" 

continuous ( HA1 4 k ")' 

in Q \ dX ep 

X 

\{HAl 4 k")'\ < 0 

at dX ep 


important” and the jump term with 5k, which could potentially compensate for the neg¬ 
ative contribution of the jump term with 5k', becomes negligible. This argument leads 
to the conclusion that the jump term with 5k' must vanish, i.e., HA1 4 k" must remain 
continuous. Since k — 0 in X e and HAl 4 ^ 0, we must have k" = 0 at dX ep . To avoid any 
doubt that this optimality condition is necessary, it is demonstrated in Appendix B that 
if the potential solutions of the localization problem are constructed with the condition 
k" = 0 at dl ep relaxed and then the minimum principle is imposed, the resulting optimum 
solution is the same as that constructed directly, with condition k" = 0 at dX ep explicitly 
imposed. 

For clarity and completeness, we compare the governing equations, boundary condi¬ 
tions, and regularity conditions for internal variable k and the two different formulations 
in Tab. 1. Conditions for the standard solution can be found in Jirasek et al. (2010) 
and Jirasek and Rolshoven (2009b). 

To simplify the following discussions, the effect of body forces will be neglected, i.e. b = 
0 in (9), which implies that the axial force 

F = EA(u' - k) (22) 


is constant along the bar. 

3. Bar With Piecewise Constant Yield Stress Distribution 

Having derived the governing equations, boundary and regularity conditions, we now 
proceed to localization analysis of a tensile test of a bar with variable initial yield stress. 
Let us consider a bar containing a weak segment of length 2 l g and an initial reference yield 
stress cr r , while the remaining parts have a larger initial yield stress ay/(1 — /3) where /3 € 
(0,1) denotes a dimensionless parameter, cf. Fig. 1. The origin of the coordinate system 
is placed into the centre of the weak segment, as we will consider symmetric solutions 
without loss of generality, i.e. X p = (— L p /2, L p /2) where L p denotes the length of the 
plastic zone X p . 
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Under these assumptions, the initial yield stress distribution is given by 


J a r for |x| < l g 

\ (Trj (1 — (3) for |x| > lg 


(23) 


with discontinuities at x — ±l g , cf. Fig. 1. The response of the bar is elastic as long as 
the stress remains below the plastic limit, and the onset of yielding occurs when the yield 
limit is attained, i.e., when F = F r with F r = Aa r = elastic limit force. 


3.1. Plastic Zone Contained in Weak Segment 

First, let us assume that the plastic zone X p is fully contained in the weak segment, 
i.e. X p C (—lg,l g ). Then, the yield condition in Eq. (12), upon substitution of A(x) = A c 
and EA{v! — n) = F = A c a c , is simplified to 

1 4 k iv (x) — k(x) = ' —- for x e Xp (24) 

H 

which is a fourth-order linear differential equation with constant coefficients and a constant 
right-hand side. It will be convenient for further analysis to convert Eq. (24) into a 
normalized form. To this purpose, we introduce 

• dimensionless spatial coordinate £ = x/l, 

• plastic strain at complete failure for the local model Kf — —a r /H , 

• normalized plastic strain K n — k/K f — —Hn/t r r , 

• dimensionless stress or load parameter (j) = F/F r , 

• dimensionless parameters describing the ratio X g = l g /l between ’’geometric” and 
material characteristic lengths, and 

• ratio X p = L p /2l between one half of plastic zone and the material characteristic 
length. 

In terms of normalized quantities, Eq. (24) is transformed into 

«n V (0 - «n(0 = 0 — 1 for £ G (-Ap, X p ) (25) 

where, for simplicity, the derivatives with respect to £ are still denoted by primes or Latin 
numerals. The general solution to Eq. (25) is 

K n(£) = Ci cos £ + C 2 sin £ + C 3 cosh £ + C 4 sinh £ + 1 — 0 (26) 


2 

'CT 0 



< 

°v/(i - /?) 

--> 


x 

Figure 1: Distribution of piecewise constant initial yield stress cro(a;) of a uniform bar. 
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where the integration constants C2 and C4 vanish due to the symmetry conditions «4(0 ) = 
0, k"'( 0) = 0. The remaining unknowns are the integration constants Ci, C 3 , and the size 
of the plastic zone X p , which are determined from the regularity conditions at the boundary 
of the plastic zone dl pj 


K n (X p ) = 0, <(A P ) = 0, <(Ap) = 0 (27) 

Substitution of the general solution (26) into the boundary conditions (27) leads to the 
set of three equations 

Ci cos X p ± C 3 cosh X p = (ft — 1 

Ci sin Xp = C 3 sinli X p (28) 

Ci cos X p = C 3 cosh X p 

Elimination of Ci and C 3 reduces the above system to a single nonlinear equation 


tan X p = tanh X p 


(29) 


Solutions of this transcendental equation can be found numerically; let us denote them 
±rrii = 0, ±3.9266, ±7.0686,... ,i G No- Of course, only positive values are of physical 
interest. For given i <G N, the integration constants can be evaluated from the Erst two 
equations of the system (28) as 


Ci = 


0-1 

2 cos rrii 


0-1 

2 cosh rrii 


(30) 


The most localized plastic strain profile is obtained for % — 1, i.e., for X p = m 1. This 
is also the standard, non-variational solution of the localization problem for a bar with 
perfectly uniform properties presented by Jirasek and Rolshoven (2009b), Eq. (40). 
Outside the plastic zone, condition (17) simplifies to 


<r c < ay or, equivalently, 0 < 1 


(31) 


Analysis of the second variation, presented in detail in Appendix A, reveals that the most 
localized solution is stable (provided that the material parameters and bar length satisfy 
a certain condition), i.e. it corresponds to a local minimum of II, while the solutions for 
i > 2 are unstable. 

So far we have assumed that the weak segment is long enough, 2 l g > 2 Imi, and thus 
the stable localized solution is not affected by stronger parts of the bar. Nevertheless, if 
the weak segment is shorter, the analysis needs to be modihed. 


3.2. Plastic Zone Extending to Strong Segments 

Let us proceed to the case when L p > 2l g , i.e. X g < mi. In this situation, Eq. (25) 
must be extended to the parts surrounding the weak segment: 


K n V (0 - «n(0 = 0-1 for -A g < £ < A g 

«n V (0 - «n(f) = 0 ± for Xg < If I < Ap 


(32) 
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yielding the general solution 


^n(0 < 


C 1 cos £ + 62 sin £ + C 3 cosh £ + C 4 sinh £ + 1 — 0 

for 

-A 9 < £ < A 

C' 5 cos £ + C e sin £ + C 7 cosh £ + C 8 sinh £ + --— — 0 

1 jJ 

for 

Xg < £ < A, 

C 9 cos £ + C10 sin £ + C u cosh £ + C V2 sinh £ + 0 

1 fJ 

for 

-X p < £ < A 


By the symmetry conditions, the integration constants C 2 and C 4 vanish again, and 

C 5 = C 9 , 

c 6 = -c 10 , 

c 7 = Cn, 

C 8 = C\2 


9 

(33) 


(34) 


The remaining unknown constants C t for i — 1,3, 5, 6, 7, 8 and the dimensionless plastic 
zone size X p can be determined from seven conditions; namely from continuity of k, 
k', HAl A n", and (HAI a k")' at £ = X g , and of k, k' , and HAl 4 n" at £ = X p . Because the 
resulting set of equations is nonlinear in X p , it is more convenient to solve the system for C* 
and 0, with X p considered as given. In other words, the loading process is considered as 
parametrized by X p instead of 0. Then, we arrive at a set of seven linear equations in the 
form 


( 

— COS Xg 

— cosh X g 

COS Xg 

sin A 


sin X g 

— sinh X g 

— sin X g 

cos A 


— COS Xg 

cosh X g 

COS Xg 

sin A 


sin X g 

sinh X g 

— sin X g 

cos A 


0 

0 

cos X p 

sin A 


0 

0 

— sin X p 

cos A 

\ 

0 

0 

cos Xp 

sin A 


cosh X g 

sinh X g 

0 ^ 


/ 

C, A 


/ 

-A_ \ 

P-1 

sinh X g 

cosh X g 

0 



C3 



0 

— cosh X g 

— sinh Xg 

0 



C 5 



0 

— sinh X g 

— cosh Xg 

0 



C 6 

= 


0 

cosh X p 

sinh Ap 

-1 



C 7 



1 

P-1 

sinh A p 

cosh Ap 

0 



Cs 



0 

— cosh A p 

— sinh Ap 

0 ) 


V 

0 / 


V 

0 / 


which can be easily solved by matrix inversion; for the sake of brevity, however, we do 
not provide the results in the explicit form. The resulting dependencies between the 
load parameter 0 and normalized plastic zone X p are depicted in Fig. 2 for several values 
of X g = mi{0.003, 0.025, 0.127, 0.255, 0.382,0.637} and for 0 = 0.5. 

Solving the system in Eq. (35) and substituting the results into Eq. (33) leads to the 
distribution of plastic strain K n . An example for X g = 0.127mi, 0 = 0.5 and the mono- 
tonically expanding plastic zone X p = mi{0.277, 0.554, 0.693, 0.776, 0.831, 0.858, 0.870} is 
shown in Fig. 3a. In Fig. 3b, the third derivative of plastic strain is depicted, satisfying 
all the regularity requirements summarized in Tab. 1. 

Integrating Eq. (33) over the length of the plastic zone provides the normalized plastic 
elongation 

u r^ p r^ p 

T ? - = «n(£)d£ = 2/ k„(£) d£ + 2 / K n (£)d£ (36) 

J Xp « 0 j Xg 
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dimensionless plastic zone size A p 


Figure 2: Piecewise constant yield stress distribution: relation between load parameter 
and plastic zone size for several values of X g , and for j3 = 0.5; for complete legend please 
refer to Fig. 4a. 


which in this simple case can be carried out analytically: 


U P _ 9 
lKf 


f3\g + Ap[(/)(1 — P) — 1] 

£-1 


+ C\ sin X g + C 3 sinh X g + C 5 (sin X p — sin X g ) 


+ C 6 (cos X g — cos A p ) + C 7 (sinh X p — sinh X g ) + C 8 (cosh X p — cosh X g ) 


(37) 


The dimensionless load-plastic elongation diagrams for fixed (3 = 0.5 and different dimen¬ 
sionless sizes of the weak segment X g are shown in Fig. 4a, and for fixed X g = 0.127mi 
with different values of [3 in Fig. 4b. These figures reflect the influence of both parameters 
on the shape of the load-displacement curve for a bar with an imperfection, where X g is 
understood as the length of the imperfection, and f3 as its magnitude. Note that regardless 
of the imperfection shape and size, the plastic response of the bar is always of the same 
slope. For shorter imperfections we observe significant hardening, while for longer imper¬ 
fections the response attains the behaviour of a perfectly uniform bar. The imperfection 
magnitude influences the load-displacement diagram in the opposite way; for large values 
of /3, the response exhibits a higher maximum force. On the other hand, in the case of 
small magnitudes the response approaches the behaviour of a uniform bar discussed in 
Section 3.1. As an alternative physical interpretation, we could consider a semi-infinite 
layer of material between two parallel planes which are mutually displaced in tangential 
direction and left unconstrained in the normal direction, so that the layer with vertically 
variable material properties is under pure shear stress. 

In order to demonstrate that we have obtained an admissible solution, we check con¬ 
dition (17), which is valid outside the plastic zone and simplifies to 

* * TTp < 38 > 

The condition can be simply verified in Fig. 2, where we have (j) < 2 for f3 = 0.5. 

Finally, in Fig. 5, we check that the energy balance (E) holds along the whole loading 
process in agrement with general results (Pham et ah, 2011, Section 2.2.3) and (Pham 
and Marigo, 2013, Property 1) for solutions sufficiently regular in time. We observe that 
the response is first elastic, i.e. the £ curve is quadratic with no dissipation Yapp, followed 
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dimensionless coordinate £ = x/l 



dimensionless coordinate £ = x/l 


(a) (b) 

Figure 3: Piecewise constant yield stress distribution: (a) evolution of plastic strain profile 
and (b) its third derivative for monotonically increasing plastic zone length X p . 



normalized plastic elongation u p /lnf 



normalized plastic elongation u p /lnf 


(a) 


(b) 


Figure 4: Piecewise constant yield stress distribution: (a) plastic part of dimensionless 
load-displacement diagram for different values of dimensionless imperfection length X g 
assuming fixed (3 = 0.5, and (b) for different values of dimensionless imperfection magni¬ 
tude /3 assuming fixed X g = 0.127mi. 
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Table 2: Physical and geometric parameters of all test examples. 


Physical parameters 

Values 

Young’s modulus, E 

60 GPa 

Softening modulus, H 

-10 GPa 

Characteristic length, l 

0.05 nr 

Reference initial yield stress, o r 

200 MPa 

Weakest cross-sectional area, A c 

0.01 nr 2 



time t 

Figure 5: Piecewise constant yield stress distribution: energy evolution for ft = 0.5; see 
Tab. 2 for the physical parameters and Eq. (39) for the loading program. 

by the evolution of the plastic strain accompanied by nonzero dissipation. For this graph, 
physical constants presented in Tab. 2 were used. Parameters that control the imperfection 
were set to X g = 0.255mi and f3 = 0.5, the total length of the bar was L = 4 Imi, and the 
evolution was parametrized by the dimensionless plastic zone length 

A p (t) = A g + t(rrii — \ g ) for t G [0,1] (39) 


4. Bar With Piecewise Constant Stress Distribution 


In all subsequent sections, contrary to the previous one, we will assume the initial yield 
stress to be constant, i.e. cr 0 (x) = oy, and will investigate the influence of the variable 
cross-sectional area resulting in spatially variable stress field distributions. As the Erst 
example, let us consider a very similar load test to that presented in Section 3, but now 
with discontinuous sectional area, i.e. a bar containing a thin segment of length 2 l g and 
sectional area A c , and with remaining thick parts of sectional area A c /( 1 — j3) where, as 
previously, [5 G (0,1) denotes a dimensionless parameter, cf. Fig. 6a. 

The corresponding stress distribution is described by 


f F/A c = cr c for |x| < l g 

\ F/[A C /(1 - 0)\ = (1 - /3)cr c for jxj > l g 


(40) 
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(a) 


(b) 



(c) 

Figure 6 : Geometries of tensile test bars corresponding to (a) discontinuous, (b) contin¬ 
uous, but not continuously differentiable, and (c) infinitely smooth stress fields. 
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cr c 
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V > 

) (1 - 0 )ct c 
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-* 
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Figure 7: Bar with piecewise constant cross-sectional area and the corresponding stress 
distribution. 


again with discontinuities at x — ±l g , cf. Fig. 7. Clearly, the case of a plastic zone 
contained in the weak segment coincides exactly with the situation presented in Sec¬ 
tion 3.1, and hence is not discussed again. Instead, let us proceed directly to the case in 
which L p > 21 g , i.e., X g < mi. 


4-1. Plastic Zone Extending to Thick Segments 

Substituting Eq. (40) into the yield condition (12) and converting the result into the 
dimensionless form, we obtain the governing equations, cf. Eq. (32): 

«n V (0 - Kn(0 = 0-1 for -\g<£<\ g 

K n V (0 - «n(f) =0-00-1 for Xg < If | < X p 

which provide the general solution 

{ Ci cos f + C2 sin f + C3 cosh f + C4 sinli f + 1 — 0 for — X g < f < X g 

C 5 cos f + Cq sin f + C7 cosh f + C 8 sinli f + 1 — 0 + 00 for X g < f < X p 

Cg cos f + C10 sin f + C'n cosh f + C'12 sinli f + 1 — 0 + 00 for — A p < f < A 9 

(42) 

By the symmetry arguments, the integration constants C 2 and C '4 are zero, and for C t , 
i = 9,..., 12, the relations (34) hold again. The remaining unknown constants Ci 

for i — 1,3, 5, 6 , 7, 8 and the dimensionless plastic zone size X p can be determined from the 

regularity conditions; for example, continuity of HA1 4 k" and {HAI a k")' at f = X g give 


(1-0K(A # -)=<(AJ) 
(i- 0 K'(a 9 -) = <(a 9 +) 
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The set of seven linear equations arising from the regularity and boundary conditions 
reads 


/ 

— COS Xg 

— cosh Xg 

cos X„ 

sin Xg 





sin X g 

— sinh Xg 

sin X„ cos X„ 





(0 — 1) COS Xg 

(1 — f3) cosh X g 

COS Xg 

sin Xg 





(1 — ( 3 ) sin Xg 

(1 — 0) sinh X g - 

sin X g cos X g 





0 

0 

cos Ap 

sin Ap 





0 

0 

sin Ap cos Ap 




V 

0 

0 

cos Ap 

sin A p 






cosh X g sinh X g 

0 \ 


f C\ 

\ 


( 

0 \ 



sinh Xg cosh X g 

0 


C 3 




0 


— 

cosh Xg — sinh X g 

0 


C 5 




0 



sinh X g — cosh X g 

0 


C 6 


= 


0 



cosh A p sinh X p 

0-1 


C 7 




-1 



sinh Ap cosh A p 

0 


C 8 




0 


— 

cosh Ap — sinh A p 

0 } 


\ ^ 

) 


V 

0 / 


and for convenience it is again solved numerically with no explicit expressions presented. 
The resulting dependencies between the load parameter 0 and the normalized plastic zone 
size A p are depicted in Fig. 8 for the same values \ g and 0 as in Section 3.2, see also Fig. 2. 

The normalized plastic strain, analytically expressed in Eq. (42), is depicted in Fig. 9a 
for A g = 0.127mi, 0 = 0.5, and for a monotonically expanding plastic zone X p = 
mi{0.277, 0.554, 0.693, 0.776, 0.831, 0.858, 0.870}. Its third derivative, presented in Fig. 9b, 
exhibits discontinuities at £ = ±A 9 and £ = ±A p . Let us note, however, that the quan¬ 
tity A(£)k" / (£) plotted in Fig. 10b has non-negative jumps only at £ = ±A P and remains 
continuous for £ e (—A p , A p ) in accordance to the discussion presented at the end of 
Section 2.1. For completeness, continuity of A1 (£)k"(£) and validity of plastic yield condi¬ 
tion (12) or plastic admissibility condition (13) can be verified in Figs. 10a and 11. 

The normalized plastic elongation, with the general expression presented in Eq. (36), 
can be again evaluated analytically: 

—L = 2 { A p [l + 0(0 - 1)] - 0/3 Xg + Cl sin X g + C 3 sinh X g + C 5 (sin A p - sin X g ) 

lK f (45) 

+ C 6 (cos Xg — cos A p ) + C 7 (sinh A p — sinh X g ) + C% (cosh A p — cosh X g )} 

Dimensionless load-plastic elongation diagrams for fixed 0 = 0.5 and for different di¬ 
mensionless sizes of the thin segment X g are presented in Fig. 12a; the influence of 0 
for fixed X g = 0.127mi with different values of (3 is shown in Fig. 12b. Notice that the 
obtained results resemble those presented in Section 3.2 for a bar with piecewise constant 
initial yield stress. However, the slope of the load-displacement diagram now strongly 
depends on X g and (3. 

Plastic admissibility condition (17) valid outside the plastic zone provides the inequal¬ 
ity already presented in Eq. (38), and can be simply verified in Fig. 8 , where for f3 = 0.5 
we require 0 < 2 . 

Finally, the energy profiles corresponding to the loading program (39) are depicted 
in Fig. 13, where we can check that the solution satisfies the energy balance (E) along 
the whole loading path. Physical constants are summarized in Tab. 2; the parameters 
reflecting the size of the thin segment were set to X g = 0.255mi and f3 = 0.5, and the 
total length of the bar was L = Alm\. 
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dimensionless plastic zone size A p 


Figure 8: Piecewise constant stress distribution: relation between load parameter and 
plastic zone size for several values of X g , and for (3 = 0.5; for complete legend please refer 
to Fig. 12a. 




dimensionless coordinate £ = x/l dimensionless coordinate £ = x/l 


(a) (b) 

Figure 9: Piecewise constant stress distribution: (a) evolution of plastic strain profile and 
(b) its third derivative for monotonically increasing plastic zone length A p . 
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0.2 


V 0 
- 0.2 


-4 -2 0 2 4 

dimensionless coordinate £ = x/l 

(a) 



dimensionless coordinate £ = x/l 

GO 



Figure 10: Piecewise constant stress distribution: (a) evolution of AL(£)k:"(£) and (b) 
A (£) K n(£) f° r monotonically increasing plastic zone length A p , and for fixed \ g = 0.127mi, 
P = 0.5. 



Figure 11: Piecewise constant stress distribution: a pictorial demonstration of plastic yield 
condition (12) valid in I p , and plastic admissibility condition (13) valid in X e for X g = 
0.127mi, \ p = 0.762mi, (3 = 0.5, and Jy = 1.5A P . 
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normalized plastic elongation u p /lnf 
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normalized plastic elongation u v /lnf 

(b) 


Figure 12: Piecewise constant stress distribution: (a) plastic part of dimensionless load- 
displacement diagram for different values of dimensionless length X g assuming fixed j3 = 
0.5, and (b) for different values of f3 assuming fixed \ g = 0.127mi. 



time t 

Figure 13: Piecewise constant stress distribution: energy evolution for f3 = 0.5; see Tab. 2 
for the physical parameters and Eq. (39) for the loading program. 
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5. Bar With Piecewise Linear Stress Distribution 


Now we proceed to a bar with a continuous but not continuously differentiable stress 
distribution, cf. Fig. 6b. The cross-sectional area corresponding to a piecewise linear 
stress distribution is specified in the form 

A(x) = , A ° l f , (46) 

lg ~ M 

where the overall length of the bar L = 2 l g now represents the supremum over all possible 
bar lengths for which the example is meaningful; note that lini.,.^ ±if/ A(x) = +oo. As in 
the previous section, A c denotes the area of the weakest cross section. Substituting A(x) 
into Eq. (12) with a 0 (x ) = ay, we obtain 


kIV{x) + 2sgnM + 


In 


\x\ 


{lg - kl) : 


;K"{x) 


— k(x) 


CT j* 

~H 



x\) for x G Z p \{0} 
(47) 


which has to be satisfied at all points inside the plastic zone with the exception of x = 0, 
where the cross-sectional area is not differentiable. At that point, we enforce continu¬ 
ity conditions of h, h', Ah" and (Ah")', recall the discussion at the end of Section 2.1 
and Tab. 1. Since A is continuous, the third condition actually reduces to continuity 
of k". After conversion to the dimensionless form, in Section 3.1, the governing equation 
transforms into 


K, 


IV 


( 0 - 


2 sgn(Q 

X 9 ~ Ifl 




'( 0 - 


(\-m 2 


«n(0-«n(0 = 0-l-0y^ for £ e (—A p , A p )\{0} (48) 


This is a fourth-order differential equation, and contrary to Eqs. (25), (32), and (41), 
it has non-constant coefficients and a non-constant right-hand side term. Although an 
analytical solution can be constructed in terms of special functions—the coefficients of the 
homogeneous part of equation (48) fulfil the so-called Calabi-Yau condition, cf. Almkvist 
et al. (2011), Eq. (3.4)—it is more convenient to solve it using the MATLAB® bvp4c 
solver, for details see Shampine et al. (2003). Again, due to symmetry conditions, it 
suffices to restrict our attention to the positive part of the plastic zone, Z+ = (0,A P ). 
Then, the boundary and symmetry conditions read 


^n(Ap) 0, K n (Ap) 0, K n (Ap) 0, and 
<(0)=0,<'(0 + ) = -<(0)/A g 


The last condition is obtained from continuity of {An")', meaning that (>l/c // ) / (0 — ) = 
(Ak")'( 0 + ), and can be derived when taking into account continuity of A, k" and skew- 
symmetry of A', k'" , i.e. A'(0 _ ) = —A'(0 + ), ^'"(0“) = —k"'( 0 + ). The solution is again 
parametrized by the length of the plastic zone A p , and for each X p the corresponding 0 is 
determined from the solution of (48) and (49). 

Before presenting the obtained results, let us briefly comment on the numerical so¬ 
lutions. Since the bvp4c solver is employed to provide only the continuous part of the 
solution on Z+, the jump condition becomes a boundary condition that is imposed directly. 
Moreover, the solver relies on a collocation method for iterative solution of boundary value 
problems with nonlinear two-point boundary conditions, and it is therefore robust with 
respect to possible changes of input data compared e.g. to shooting-based strategies used 
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Figure 14: Piecewise linear stress distribution: (a) plastic part of load-displacement 
diagram, and (b) relation between load parameter and plastic zone for X g = 
{1.019, 2.037,4.075, 8.150}mi. The maximum value of the load parameter 0 decreases 
with increasing X g . 


in Jirasek and Zeman (2015). For further details, we refer to Shampine et al. (2003), 
Section 3. 

The plastic part of the load-displacement diagram is depicted in Fig. 14a, where 
the dimensionless load parameter 0 is plotted against the dimensionless plastic elon¬ 
gation Up/lKf. The initial part of the diagram is vertical, as in the previous examples, 
since only the clastic deformation evolves for F < F r . At the onset of yielding, the load 
parameter first steeply increases and only later decreases. Complete failure is attained at 
larger elongation in comparison with the non-variational formulation analyzed in Jirasek 
et al. (2010). Several values X g = {1.019, 2.037,4.075, 8.150}mi are reported (from top to 
bottom), to reflect the effect of spatial variation of the sectional area; lower values of X g 
correspond to a stronger variation of the sectional area and lead to higher peak loads. 
Note that for the standard formulation, the overall elongation at failure is always the 
same, while for the variational approach it depends on X g . 

Fig. 14b captures the evolution of the plastic zone size. The load parameter 0 is plotted 
against the dimensionless plastic zone size X p , obtained again for several values of X g . We 
can infer from the figure that the plastic zone evolves continuously and monotonically 
from the weakest section to its full size. The length at complete failure is almost the same 
as for the standard solution, plotted by dashed curves. 

The evolution of the plastic strain profile and of its third spatial derivative is depicted 
in Fig. 15 for X g = 1.019mi and several values of X p = mi{0.255, 0.637, 0.764, 0.891, 0.998}. 
First, during the early stages of plastic evolution, the standard and variational solutions 
are almost the same, but at later stages, the differences grow significantly. Contrary to the 
standard solution, k'" is discontinuous for the variational solution at £ = 0, where [[/v /// ] 0 = 

-2k"(0 )/l„. 

In the elastic zone Z e , condition (17) for an admissible solution reduces to 

0 < —^ (50) 

1 
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dimensionless coordinate £ = x/l 


(a) 


(b) 


Figure 15: Piecewise linear stress distribution: (a) evolution of plastic strain profile and 
(b) third derivative of plastic strain for monotonically increasing plastic zone length X p . 




time t time t 


(a) standard 


(b) variational 


Figure 16: Piecewise linear stress distribution: energy evolutions; see Tab. 2 for the 
physical parameters and Eq. (52) for the loading program. 
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where it is sufficient to verify £ = X p . For the data in Fig. 14b we obtain 


14 

**I3T = 3 < 51 > 

and the condition is satisfied. For X p —> X g , however, the right-hand side in (50) converges 
to +oo showing that plastihcation of points close to physical boundary d£l would require 
a very strong growth of (ft. 

Figure 16 depicts the energy balance (E) for the standard solution (Fig. 16a) and for 
the energetic solution (Fig. 16b) obtained for the data presented in Tab. 2. Further, we 
have used X g = 1.02mi and have parametrized the localization process through 

X p (t) = m\t for t G [0,1] (52) 

It is worth noting that, for the standard solution, the work done by external forces JqV(s) ds 
is out of balance with £ + Varx>, while the variational approach delivers an energy- 
conserving process. 


6. Bar With Quadratic Stress Distribution 

As the final example, we shall present the most regular case of quadratic stress distri¬ 
bution possessing continuous derivatives of an arbitrary order, see Fig. 6c. The function 
describing the cross-sectional area then reads 

A c l 2 a 

A(x) = (53) 

Lg ~ X Z 


where the overall length of the bar L = 2 l g is, in analogy to Eq. (46), understood as the 
supremum over all possible bar lengths for which the example is meaningful; note that 
again, lim x _^-t/ A(x) = +oo. Upon substitution into the yield condition (12) with a 0 (x) = 
ay, we get the governing equation 


4x 


2 (P + 3x 2 


l 4 K W (x) + r) 2 k"'(x) + ^ K"(x) 


1 2 9~x‘ 


ill ~ 


- «(®) = ^ _ ^ ( 1 _ 7 ^ ) fori6l p 


(54) 


which can be converted into the dimensionless form 


jv ft \ i J"(t\ , 2 ( A s + 3 ^) 


C(0 + 




+ 


(A 2 - £ 2 ) 2 


K n(0 - Kn(0 =(ft-l- for £ G (-Ap, X p ) (55) 


As in the previous case, we will employ a numerical solver, since the governing equation 
is even more complicated. Owing to symmetry requirements, the solution will again be 
constructed in the positive half of the plastic zone Z+, with boundary and symmetry 
conditions 

K n( Ap) = o, K n (x p ) = 0, K n (Xp) = 0 

<(0) = 0,<(0) = 0 1 J 

In contrast to (49), the last condition is now simpler since A 1 is continuous. 

The solution has been computed for several values of X p . In Fig. 17a we notice 
that the plastic zone evolves continuously and monotonically; particular plastic strain 
profiles together with their third derivatives are depicted in Fig. 18 for X g = 1.273mi 
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Figure 17: Quadratic stress distribution: (a) relation between load parameter 
and plastic zone size, (b) plastic part of load-displacement diagram for X g = 
{1.019,1.273, 2.547,12.734}mi. The maximum value of the load parameter 0 decreases 
with increasing X g . 


and X p = mi{0.255, 0.636, 0.764, 0.891, 0.968, 0.998}. Comparing the results presented in 
Fig. 15 with the results in Fig. 18, we notice that for the quadratic stress distribution the 
differences between the standard and variational solutions are less pronounced. Due to a 
higher smoothness of the solution, the load-displacement diagram presented in Fig. 17b is 
almost linear, only with a slight hardening followed by the softening branch. Differences 
between plastic displacements at failure, i.e. for 0 = 0 , are also somewhat less distinct. 

Hardening effects for the variational formulation are systematically stronger in com¬ 
parison to the standard formulation; moreover, for small values of X g , the differences 
are more obvious. This effect has already been explained in Section 2.1, see Eq. (12), 
Tab. 1 and the discussion therein. Recall, nevertheless, that the two formulations differ in 
two terms with higher-order derivatives of the sectional area, neglected for the standard 
formulation. For decreasing magnitudes of A!{x) and A!'{x), the variational formulation 
approaches the standard one; the limit case is presented in Section 3.1, where the two 
solutions coincide. Let us note that for an infinitely differentiable exponential stress dis¬ 
tribution, considered in Jirasek et al. (2013), we would also obtain significant differences 
for X g small enough. 

Substituting expression (53) into inequality (17) leads to 



which should hold inside Z e , i.e. for all |£| > X p . A closer inspection of Fig. 17a reveals 
that the condition is satisfied. 

Energy balances for the standard and variational formulations are shown in Fig. 19 
using the data from Tab. 2 and loading program in Eq. (52). The total length of the 
bar was 2X g with X g = 1.019mi. Again, for the standard formulation we notice a slight 
violation of condition (E). 
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dimensionless coordinate £ = x/l dimensionless coordinate £ = x/l 

(a) (b) 

Figure 18: Quadratic stress distribution: (a) evolution of plastic strain profile and (b) 
third derivative of plastic strain for increasing X p . 




time t 


time t 


(a) standard 


(b) variational 


Figure 19: Quadratic stress distribution: energy evolutions; see Tab. 2 for the physical 
parameters and Eq. (52) for the loading program. 
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7. Summary and Conclusions 

We have presented one-dimensional localization analysis of a softening plasticity model 
regularized by a variational formulation with a fourth-order gradient enrichment. The 
main results are summarized as follows: 

1. Using a consistent variational approach, we have derived the description of a one¬ 
dimensional gradient plasticity model which provides not only the appropriate dif¬ 
ferential equation, representing the yield condition inside the plastic zone, but also 
appropriate forms of boundary and jump conditions at the elasto-plastic interface. 

2. On the basis of the derived conditions that follow from the variational principle, two 
problems with discontinuous data (a bar or layer with discontinuous yield stress and 
a bar with discontinuous cross-sectional area) have been investigated. These exam¬ 
ples are amenable to analytical solution and have provided physically reasonable 
results. 

3. Two additional examples have been analysed, one with continuous but not contin¬ 
uously differentiable data and the other with smooth data. Numerical solutions 
have been constructed and compared to the alternative non-variational formula¬ 
tion, demonstrating that the variationally consistent formulation leads to higher 
peak loads and elongations at structural failure. 

4. We have also investigated the influence of various data on the evolution of the 
plastic zone and on load-displacement diagrams for all four prototype problems. 
It has been shown that the plastic zone monotonically expands from the weakest 
section of the bar. In spite of the softening character of the material model, the 
structural response exhibits first hardening after the onset of yielding, later followed 
by softening. Such a behaviour is related to a gradual expansion of the plastic zone 
to stronger segments of the bar. 

5. Further, it has been demonstrated that the solution corresponding to the variational 
formulation satisfies the energy balance along the evolution path of the localization 
process. Contrary to that, the standard formulation exhibits systematic lack of 
balance in the sense that the sum of elastic and dissipated energies exceeds the 
work done by external forces. In all cases, however, the dissipated energy remains 
finite and non-zero. 

6. The variational approach is based on the non-negative first variation of the func¬ 
tional. Solutions corresponding to local minima of the functional have to satisfy 
this condition and moreover have to be stable. Hence, an analysis of the second 
variation providing some explicit requirements on physical constants is presented 
in Appendix A for the simplest case of a bar with perfectly uniform data. 

Let us note that although we have not employed variable elastic or plastic moduli, the 
variational approach is perfectly suited to handle problems with discontinuities in such 
data. 

The presented framework can also be extended to higher dimensions using e.g. the 
von Mises yield function with isotropic softening. Then, the free boundary conditions for 
the scalar cumulative plastic strain at the elastoplastic interface are analogous to those 
for k for sufficiently smooth fields only. However, in the multi-dimensional setting the 
required regularity is difficult to establish since, for instance, the embeddings W 1,2 C C° 
or W 2 ' 2 C C 1 no longer hold; its rigorous investigation is beyond the scope of the present 
work. 
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Appendix A. Second Variation and Stability Conditions 

The analysis presented in the main part of the paper has been based on the condition 
of non-negative first variation of the energy functional fl. This condition is, however, 
only a necessary one, yet is not sufficient to ensure that the solution is a local minimum, 
thus that it is energetically stable. Therefore, in this section we discuss the behaviour 
of the energy functional in the vicinity of a solution (u,k) satisfying (11)—(13) , drawing 
inspiration from a related analysis by Jirasek et al. (2013). 

In particular, we will investigate the second variation of the energy functional given by 
Eq. (7). Our objective is to show that the second variation is positive for all those nonzero 
admissible variations Su and 5 k for which the first variation <511 vanishes, to ensure that 
the solution (u, k) is stable. Overall procedure will be described in six steps for better 
clarity. In Appendix A.l, we start with the elimination of the displacement field from the 
functional II in order to simplify stability conditions discussed in Appendix A. 2, where 
the corresponding eigenvalue problem will be derived. In Appendix A.3 and Appendix 
A.4, we discuss even and odd eigenfunctions and specify the requirements on physical 
and geometric parameters leading to stable responses. A discussion of larger plastic zones 
for a uniform bar is presented in Appendix A. 5. In Appendix A.6, we summarize our 
developments and briefly comment on problems with non-uniform cross-section area. 


Appendix A.l. Condensation of Displacement Field 

In the first step, we can simplify the problem by eliminating the displacement field and 
constructing a reduced functional that depends on the plastic strain field only. Indeed, the 
original functional (7) can first be minimized with respect to the displacements at fixed 
plastic strains. Integrating the already derived optimality condition (11) and assuming 
vanishing body forces 6, we obtain 

u'{x) = f + k(x) (A.l) 

EA(x) 

where F is the integration constant, physically corresponding to the axial force, cf. 
Eq. (22). Integrating again and taking into account the boundary conditions implied 
by (2a), we can express the normal force as 


F — K e I u — / k(x) da: 


where u(t) = UD(dO,R,t) — UD(dCLi,,t) is the prescribed total elongation, and 

1 


Ke = 


f dx 

In BA(x) 


(A.2) 


(A.3) 


is the elastic stiffness of the bar (reciprocal value of the elastic compliance). Based 
on (A.1)-(A.2) and on the assumption of no body forces (6 = 0), the original functional (7) 
can be reduced to 

H*(k) — -^(u — f /fda:^ + - f HA(k 2 — I^k" 2 ) dx + f Aafkdx (A.4) 

2 \ Jn J 2 J n J n 


Note that the first term represents the elastically stored energy, expressed as a functional 
dependent on the plastic strain only. The expression in the parentheses is the elastic 
elongation, written as the difference between the total elongation u and the plastic part 
of elongation u p which was introduced in (36). 
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Appendix A.2. Second Variation 

The first variation of the reduced functional fl* is given by 


5U*(k; 5k) =—K e [ u — I ndx ) / 5Kdx + / F[ A(k5k — I a k" 5k") da; + / Aa 0 Sn dx 

'n Jn Jn 

(A ' 5) 

which exactly corresponds to (8) with v! and 5u' expressed according to (A.1)-(A.2) and b 
set to zero. The second variation of the reduced functional TT is given by 


5 2 L[*(5k;) = K e ( dndx) + / HA(5k 2 - F5k" 2 ) da; 


(A.6) 


and is independent of k because fl* is quadratic. 

We would now like to prove that the quadratic form 5 2 II* is positive definite in the 
space of all those admissible variations 5k for which 511* (/?;£/?) = 0. Note that k is 
not arbitrary but represents the solution of the problem described by conditions (14)— 
(16), which were derived from the requirement 5II(«, k\ 5u, 5k) > 0 for all admissible 
variations 5u and 5k (and which could alternatively be derived, in a slightly different 
format, from the condition 511* (/ c ; 5k) > 0 for all admissible variations 5k). Integrating 
the term with 5k" in (A.5) twice by parts and exploiting the properties of function k (in 
particular, the fact that k = 0 in X e ), we can rewrite the first variation of II* as 


5II*(k; 5k) 



HAk - (. HAI a k ")" + Aa 0 - K e 




5 k dx 


+ 



Aa 0 - K e 



5k dx 


- HAl A k" n5 k' + Y( HAl * K ")' nSK 

dl ep dl ep 


(A.7) 


where n denotes the outer normal to X p . Taking into account that k satisfies the yield 
condition (12) in which the right-hand side corresponds to the normal force F , which 
is in turn equal to the expression given in (A.2), we can show that the expression in 
square brackets in the first integral in (A.7) vanishes. The first sum in (A.7) also vanishes 
because k" = 0 on dX ep . Consequently, (A.7) reduces to 


5n*(/c; 5k) 



5 k dx + Y( HAl4 K")'n5K 

dX ep 


(A.8) 


Let us also recall that (i) the variation 5k is nonnegative in X e UdX ep ] (ii) the expression 
in the square brackets in (A.8) is nonnegative in Z e , see condition (17); and (iii) the 
product (HAl 4 K")'n is nonnegative on dX ep , see the last line in Table 1. To proceed 
further, we need to assume that the expression in the square brackets in (A.8) is strictly 
positive (not just nonnegative) in X e . This assumption means that the normal force in the 
entire elastic domain is below the yield limit, which is true for all the localized solutions 
presented in this paper. Consequently, for those admissible variations 5k that are not 
identically zero in X e , the integral in (A.8) is positive and thus the first variation 511* is 
positive, because the second term on the right-hand side of (A.8) is always nonnegative. 
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Based on the foregoing analysis of the first variation, we can see that positive definite¬ 
ness of the second variation h 2 fl* needs to be proven for those variations 8k that vanish 
in X e , and the integration domains in (A.6) can be restricted to X p . To preserve continuous 
differentiability of 8k in hi, its value and first derivative on dX ep = dX v must be zero. To 
simplify notation, we will use from now on v instead of 8k. The stability condition is then 
satisfied if there exists some a > 0 such that 

K e (f vdx ) + f HA(v 2 — l A v" 2 )dx > a||u|| 2 for all v e V (A.9) 

\Ji P J Jip 

where ^ 

V = W r f' 2 (Z p ) = [v e W V \Z P ) I v\ dIp = 0, t-V, = 0} (A. 10) 

and where ||u|| denotes standard L 2 (Z p )-norm of a function v. 

Since K e > 0 and HA < 0, it is instructive to rewrite inequality (A.9) as 

Ii e ( f vdx\ — f HAl A v" 2 dx>— f HAv 2 dx + a\\v\\ 2 (A.11) 

y JJ JJ 

or, in a more abstract form, as 

(v,v) A > (v,v) B + oi\\v \\ 2 (A.12) 

where (u,i’)A and ( v , v)b are quadratic forms corresponding to symmetric bilinear forms 


(v,w)a = K e vdx wdx— HArv"w"dx 


(' v,w) B = 


HAvw da; 


(A. 13) 
(A. 14) 


Since form ( v,w)b is positive definite and form (v, w)a is at least positive semidehnite (in 
fact it is positive definite), condition (A.12) can be reformulated as 


X _ . f (■ v,v)a - n 

-^min — inf / \ ^ 1 

nev\{ 0 } {v,v)b 

Indeed, if (A. 15) is satisfied, then we have (for all v G V) 


(A.15) 


(■ V,v)a > A min (v,v) B = ( V,v)b + (A mi „ - l)(v,v) B > ( V,v) B + (A m in ~ 1 )(~HA) 1 


(A.16) 

and thus inequality (A.9) is satisfied with a = (A min — 1 )(—HA) min , where (—HA ) min > 0 
is a lower bound on —HA over X p . 

The fraction in (A.15) is the well-known Rayleigh quotient, and A m i n is the minimum 
eigenvalue obtained from the eigenvalue problem 

(n, w)a — A(u, w)b — 0 for all w G V (A.17) 

Note that the rigorous treatment of (A.15), its transition to (A. 17), and relation to the 
Poincare inequality can be found in Sebestova and Vejchodsky (2014). 

To finish the stability analysis, it is necessary to find the eigenvalues A for which (A. 17) 
has a nontrivial solution v and check that they are greater than 1. In general, this would 
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need to be done numerically. However, it is instructive to construct an analytical solution 
for the simplest case of a uniform bar, with EA = const., HA = const., and K e = EA/L. 
The eigenvalue problem (A. 17) then reads 


EA 

~T~ 


vdx wdx — HAI a / v"w" dx + A HA / vw dx 


0 for all w £V (A. 18) 


Integrating by parts and exploiting the boundary conditions w = 0 and w' — 0 on dl p , we 
can construct the corresponding eigenvalue problem in terms of differential (and in this 
case also integral) operators: 


- vdx- Hl A v IV + A Hv = 0 (A. 19) 

L Ji P 

This is an integro-differential equation, which can be for convenience converted into the 
set of two equations, 

74 IV 4 /~i 

l v — uj v = Co 
E f 

HL i t Vix = a 

Here, cu = A 1 / 4 , just to simplify the subsequent derivations, and Co is an auxiliary un¬ 
known. 

The advantage of the transformation of (A. 19) into (A.20)-(A.21) is that (A.20) is 
a standard linear differential equation with constant coefficients and constant right-hand 
side, and its general solution is easily expressed as 

/ N Co „ LUX . UX _ , UJX _ . , LUX . . . 

v{x) =-- + Ci cos ——h C 2 sm ——h C 3 cosh ——h C 4 smh — (A.22) 

UJ L l L L 

Recall that the plastic zone for the most localized solution in a uniform bar corresponds to 
the interval X p = (—mil, mil) where m\ = 3.9266 is the smallest positive solution of the 
equation tanx = tanhx. The solution v £ V must satisfy boundary conditions u(±mR) = 
0 and i/(±mil) = 0, which provide four equations for five unknown constants, Co to C 4 . 
The fifth equation is obtained from (A.21). By simple row manipulations, the resulting set 
of five homogeneous linear equations can be decoupled into two independent sets, which 
read 


(A. 20 ) 
(A.21) 


and 


— oj 4 Co + Ci coscami + C 3 cosh cum \ = 0 

—Ci sin iumi + C 3 sinh umi = 0 

Co + Ci sin iumi + C 3 sinh umi = 0 


mi HLou A 
+ 2 El ) 


(A.23) 
(A.24) 

(A.25) 


C 2 sin cumi + C 4 sinh cami = 0 
C 2 cos iumi + C 4 cosh = 0 


(A.26) 
(A.27) 
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Figure A.20: Relation between dimensionless parameter —HL/El and eigenvalues of the 
first kind. 


Appendix A.3. Even Eigenfunctions 

Let us first examine (A.23)-(A.25), having a nontrivial solution if and only if 
-4 ( m i HLu\ 

2u sinwmismhwmi — — + (cos wmi smh urnii + sinwmi coshwmi) = 0 

\ oj 6 2 El J 

(A.28) 

which is equivalent to 


or 


mioj 


2 sin umi sinh urnii 


cos urn \ sinh urn \ + sinumii cosh com i 


HL 

w 


(A.29) 


Positive solutions of this transcendental equation, ny , z = 1,2,..., correspond to eigen- 

values A^ = ) , i = 1, 2 ,.... Superscript ^ means that we are dealing with solutions 

of the “first kind”, resulting from singularity of (A.23)-(A.25), with nonzero constants Co, 
Ci and C 3 and with C 2 = C 4 = 0. Consequently, the corresponding eigenfunctions are 
even. The eigenvalues cannot be expressed analytically, but one can characterize them 
graphically, by plotting the left-hand side of (A.29) as a function of A = ca 4 , as shown in 
Fig. A.20. The vertical axis corresponds to the dimensionless ratio —HL/El. For a given 
bar, this ratio is fixed and the intersections of the corresponding horizontal line with the 
graph provide the eigenvalues for the considered case. 

From Fig. A.20 it is clear that if —HL/El is at or above a certain critical value, 
the smallest eigenvalue of the first kind is smaller than or equal to 1 and the stability 
condition is violated. The critical value is obtained simply by evaluating the left-hand 
side of (A.29) for uj — 1, and turns out to be equal to 2 (mi — tan mi) = 5.8548. Thus 
the resulting stability condition reads 


HL 

Jn 


< 2 (mi 


tan mi) 


(A.30) 


This is a constraint that involves the ratio between the absolute value of the softening 
modulus and the elastic modulus, —H/E, and the ratio between the bar length and the 
material characteristic length, L/l. Stability (under displacement control, i.e., with u 
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dimensionless coordinate £ = x/l 



dimensionless coordinate £ = x/l 


(a) (b) 

Figure A.21: Eigenfunctions corresponding to the smallest eigenvalue (a) of the first kind, 
(b) of the second kind. 


prescribed) is compromised for long bars (large L/l) made of materials with pronounced 
softening (large —H/E). The eigenfunction 

x x 

cos — cosh — 

2 cos rri\ 2 cosh mi ’ 

corresponding to the smallest eigenvalue of the first kind in the critical case when = 1 
(i.e., when — HL/El = 2 (mi —tanmi)) is plotted in Fig. A.21a. Up to an arbitrary scalar 
multiplication factor, it is identical with the actual localized solution k(x) described in 
normalized form by equations (25) and (30). The critical case corresponds to a snapback 
point of the load-displacement diagram, i.e., to a point at which the tangent becomes 
vertical and the amplitude of the plastic strain profile can grow at constant total elongation 
of the bar while the equilibrium equation and the yield condition remain satisfied. 


x G I p = (—mil, mil) (A.31) 



Appendix A.4- Odd Eigenfunctions 

We still need to examine another set of eigenvalues, referred to as eigenvalues of the 
“second kind”. Equations (A.26)-(A.27) have a nontrivial solution if and only if 


sin cam i cosh cam i — coscami sinhcami = 0 


which is satisfied for 


(2) mi • i o 

u>i = —, i = 1 , 2 ,... 

m i 


The corresponding eigenvalues and eigenfunctions of problem (A. 19) are 

4 

k ( 2 ) 


( 2 ) 

V> ' = 


(cafV = 

fm 

V 1 J 

\mi 

mix 

sinh 

sin-- 

mil 



miX 

mil 


sm m; 


sinh m,- 


, x e = (— mil, mil) 


(A.32) 
(A.33) 


(A.34) 
(A.35) 


( 2 ) 

The smallest eigenvalue of the second kind, A; , is equal to 1, while all the others are 

( 2 ) 

greater than 1. All eigenfunctions of the second kind are odd, and the eigenfunction v\ ' 
corresponding to the smallest eigenvalue is plotted in Fig. A.21b. 
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dimensionless coordinate £ = x/l dimensionless coordinate £ = x/l 

(a) (b) 

Figure A.22: Eigenfunctions corresponding to the smallest eigenvalue (a) of the first kind, 
(b) of the second kind. 


Appendix A.5. Note on Larger Plastic Zones 

The analysis presented so far for a uniform bar has referred exclusively to the localized 
solutions with the smallest possible plastic zone of size L p = 2mil. The original set of 
conditions (14)—(16) , from which the plastic strain distribution k(x) was determined, may 
admit other solutions, with larger plastic zones. Indeed, the plastic zone size follows from 
equation (29), which has positive solutions X p = m *, z = 1,2,..., corresponding to plastic 
zone sizes L p = 2 md, i = 1,2,.... In Section 3.1 we decided to focus on the most 
localized solution with i — 1. If a uniform bar is long enough, there exist other solutions, 
corresponding to i > 2 . An example for i = 2 is plotted in Fig. A.22a. Such solutions 
satisfy conditions that follow from non-negativeness of the first variation, but they violate 
the stability condition and thus cannot physically occur. To see that, consider the case 
of i — 2, i.e., X p = (—m 2 Z,m 2 Z) where m 2 = 7.0686. Condition (A.32) is then replaced by 

sin u;m 2 cosh u;m 2 — cos wm 2 skill cum 2 = 0 (A.36) 


Eigenvalues of the second kind are given by 



(A.37) 


and the minimum eigenvalue, = (mi/m 2 ) 4 , is smaller than 1. This indicates that 
condition (A.15) is violated and stability is lost. Fig. A.22b shows the eigenfunction 
corresponding to the minimum eigenvalue. 


Appendix A.6. Summary 

In summary, for a uniform bar (or sufficiently long uniform weakest segment of a 
bar) we have found that eigenvalues of the first kind are greater than 1 if the material 
constants and bar length satisfy a certain constraint, but the smallest eigenvalue of the 
second kind is always equal to 1. Strictly speaking, this would mean that stability cannot 
be guaranteed, because inequality (A.9) cannot be satisfied with any positive a but only 
with a = 0. The second variation is then non-negative but not positive definite. This is 
related to the special nature of the localization problem for a uniform bar, which actually 
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/ 9 \ 

Figure A.23: Eigenvalue of the second kind, A), as a function of a G [0, for cross- 
sectional areas corresponding to piecewise linear and quadratic stress distributions, see 
Eq. (A.38). 


admits infinitely many localized solutions characterized by the same energy level. Indeed, 
the precise position of the localized zone in a perfectly uniform bar remains undetermined, 
and the considered solution (centered in the middle of the bar) can be arbitrarily shifted 
without modifying the total energy. In a real bar, the position of the localized plastic 
zone would be affected by imperfections of the geometry and material data (sectional 
area, initial yield stress, etc.) For an ideal, perfectly uniform bar, the localized plastic 
zone can form anywhere and the corresponding solutions are neutrally stable. It is not 
by chance that the eigenfunction v ^ corresponding to eigenvalue A^ = 1 is actually a 
scalar multiple of the spatial derivative of function k which describes the plastic strain 
distribution. Adding an infinitely small multiple of k' to k corresponds to an infinitesimal 
horizontal shift of the localized plastic strain profile. Note that adding a finite multiple 
of k' would result into violation of the constraint k > 0 near one boundary of the plastic 
zone and is thus inadmissible. Shifted plastic strain profiles form a parametric family 
of functions at the same energy level which does not correspond to a linear manifold in 
the functional space. Convex linear combinations of two selected members of this family 
correspond to higher energy levels and if the difference of two such functions is considered 
as a variation 5k, the corresponding first variation 511* is positive. This is in agreement 
with our conclusion that variations which do not vanish everywhere in the elastic zone 
lead to positive 511* and thus do not need to be considered in the stability analysis based 
on the second variation. 

Before closing Appendix A, let us demonstrate the presented approach for a bar with 
a variable cross-sectional area. In particular, we verify A) > 1 by numerically solving 
the eigenvalue problem (A. 17) for a family of cross-sectional areas 


Ai(x) 

A q {;x) 


1 + 
1 + 


a\x\ 


1 — a\ 

\x\ 

(ax 

) 2 


1 — (ax) 2 


for x G l p = (—mil, mil) 


(A.38) 


where a G [0, ^A), cf. also Eqs. (46) and (53), corresponding to piecewise linear and 
quadratic stress distributions. The obtained results are depicted in Fig. A.23 where we 
verify that for a —> 0 we get A) ; —> 1, and that A) 7 > 1 for a ^ 0, as expected. 
This numerical result confirms that nonuniformity of the cross-section area A(x) has a 
stabilizing effect. 
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Appendix B. Stability of the Homogenenous Boundary Condition k " — 0 


In order to verify that [//A/ 4 ^'] = 0 is the correct optimality condition at the clasto- 
plastic interface, we construct a family of solutions without taking this condition into 
account, and then prove that the energy-minizing state coincides with the solutions pre¬ 
sented above. 

For simplicity, we consider a uniform bar described by (25) with a symmetric general 
solution presented in (26), where C 2 = 0 and C 4 = 0 , and where the boundary conditions 
read 

«n(A p ) = 0, <(A„) = 0 (B.l) 

We have two conditions, but three unknowns, Ci, C 3 and A p . The last condition in (27), 
i.e. k"(\ p ) = 0, is not imposed, but will be justified by direct energy minimization. 
Substituting (26) into (B.l) for X p > 0 yields the solution 


K, 


i(f) = !- </>- 


cos Ap sinh X p + cosh A p sin A p 


(sinh Ap cos f + sin X p cosh £) (B.2) 


In order to make our exposition more readable, the subsequent derivations are struc¬ 
tured into four steps. In Appendix B.l, we first normalize the functional II* from Eq. (A.4) 
in order to reparametrize it in terms of the dimensionless solution (B.2). Such a formu¬ 
lation makes it relatively easy to demonstrate that the homogeneous interface conditions 
correspond to a saddle point of the reduced energy functional, both under fixed axial force 
or prescribed displacements described in Appendix B.2 and Appendix B.3. In Appendix 
B.4 we finally demonstrate that the saddle points are energy minima. 


Appendix B.l. Normalization of the Reduced Functional 

We search for the minimizers of the energy functional II*, which after normalization 
takes the form 


n(fc„) = 0 




— K n + K 'n + 2 Kn) d£ 


(B.3) 


where we have denoted 


H2lg 

In (B.3), u = u/uo denotes the dimensionless elongation, u 0 = 2l g a r /E = 2 l g e 0 , 
Al g of/E denotes the reference energy. Since F = Aa r <f> , we rewrite (A.l) as 


(B.4) 


and 0 = 


v! — K = 


F 

~EA 



(B.5) 


Integration over the plastic zone X p and conversion to the normalized form provides 


= u — 6 


K n d£ 


(B- 6 ) 


Ar? 


see also (A. 2 ), which can be introduced into the energy functional (B.3) and furnishes us 
with the relation 

n(re„) = 00 2 + Ip9 f (-/4 + «n 2 + 2 «n) d£ (B.7) 

J Ap 

Two situations can now be investigated: minimization under fixed axial force 0, or under 
prescribed displacement u. Note that, for a prismatic bar with uniform properties in 
inelastic regime, we have 0 < 1, which can be verified in Figs. 4b and 12b for f3 —» 0 or 
Eq. (31), and that u > 1 directly from its definition. 
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Appendix B.2. The case of fixed <f> 


Direct differentiation and integration of (B.2) provides 


[ K n d£ = 2(1 — 0)(Ap - 2 a) 

J A p 

(B. 8 a) 

(~ K n + K 'n) d£ = 2(1 — 0 2 )(A p — 2a) 

J Ap 

(B. 8 b) 

where we have denoted 


. ^ ^ sinhApSinAp 

cos Ap sinh X p + cosh X p sin X p 

(B.9) 

Substituting from (B. 8 ) into (B.7) gives after some algebra 


S(Ap) = 0 {0 2 + 20(1 - 0 2 )[Ap - 2a(Ap)]} 

(B.10) 


from which we obtain the stationarity condition (the prime now denotes the derivative 
with respect to A p ) 


n'(Ap) = 200(1 - 0 2 ) —[A p - 2a(A p )] = 0 


(B.ll) 


which reduces to (29) and hence X p — rri\. Taking the second derivative provides 
d 2 


dA 2 


[A p — 2a(A p )] 


\p=m\ 


d 2 

= —2 ^2 a (Ap) 


= 0 


(B.12) 


\ p =mi 


From Fig. B.24 and Eq. (B.12) we deduce that X p = mi is a saddle point of IT. 


Appendix B.3. The case of fixed u 

Introducing (B. 8 a) into Eq. (B. 6 ) provides 


u — 20[A P — 2ct(Ap)] u — 1 

1 - 26*[Ap - 2a(A p )] = " D{X p ) 


(B.13) 


where we have denoted D(X p ) = 1 — 29[X P — 2a(X p )\. Substituting (B.13) 
gives 


n(Ap) = 0 


[1- D(X p )}(l-2u) + u 2 

D(\) 


into 


(B.10) then 
(B.14) 


which is the normalized energy functional under prescribed fixed dimensionless elonga¬ 
tion u. Minimization with respect to X p yields 


n'(Ap) 


- 0(1 — u ) 2 d 
D 2 (X P ) IX P 


D( Ap) = 200 


(1 — u ) 2 d 
~DfiX p )~dX p 


[Ap 


2a(Ap)] = 0 


(B.15) 


The condition ^-[A p — a(A p )] = 0 reduces again to Eq. (29), cf also (B.ll). The second 

derivative of II provides again condition (B.12) showing that the solution is also at a 
saddle point. 
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Figure B.24: n(A p )/n(mi) for fixed 0 = 0.5, and n(A p )/n(mi) for fixed u = 3 as functions 
of A p G [ni,mi] and for a bar with uniform properties. 


Appendix B.4- Energy minima 

To verify that X p = mi is actually the minimum, we recall the geometric constraint k"(A p ) > 

0 resulting from the requirements that n(X p ) = Q,k'(X p ) = 0, ft(£) > 0 , k(£) = 0 

for £ ^ (—A p , A p ), and a Taylor series expansion in X p near the boundary point A p . 
Consequently, the constraint ac"(A p ) > 0 gives for the general solution in (B.2) condi¬ 
tion sinhA p cosA p — sin A p cosh A p > 0, which provides the constraint A p G [ni,mi\. For 
the definition of mi see Eq. (29) and the discussion below; analogously we define n t as 
solutions of tan rij = — tanhrij, which leads to n ±i = 0, ±2.3650, ±5.4978,..., i G N 0 . 
Since 

n'(Ap) = n'(Ap) = Cn4~[ X r ~ MK)} 

dAp 

(coshA p sinA p — cosA p sinhAp ) 2 f < 0 for \ p G [ni,mi) 

n (cosh Ap sin X p + cos X p sinh A p ) 2 1 = 0 for mi 

(B.16) 

where C'u > 0 is a A p -independent constant, cf. Eqs. (B.ll) and (B.15), we conclude 
that the minimum is attained for X p = mi. This finding can be visually verified in 
Fig. B.24 constructed for the data from Tab. 2, for L = 2l g = 4Imi, and for the load 
parameters 0 = O.5oru = 3. 

In conclusion, for both cases, i.e. either for fixed 0 or u, the boundary condition /c"(A p ) = 

0 is indeed the optimal one and the solution is located at a saddle point, which is at the 
same time the boundary of the admissible set [ni,mi]. 
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